classdef PhysicalVal < handle
    properties
        mean_val
        stderr
        symbolic
        units
        derived_from
    end
    
    methods
        function obj = PhysicalVal(varargin)
            if nargin<1
                obj = PhysicalVal(0,0,1);
            elseif nargin==3
                val = varargin{1};
                err = varargin{2};
                unit = varargin{3};
                allowed_symbols = 'a':'z';
                str_length = 12;
                name = allowed_symbols(randi(numel(allowed_symbols), [1 str_length]));
                prev_vars = containers.Map;
            elseif nargin==4
                name = varargin{1};
                val = varargin{2};
                err = varargin{3};
                unit = varargin{4};
                prev_vars = containers.Map;
            elseif nargin==5
                name = varargin{1};
                val = varargin{2};
                err = varargin{3};
                unit = varargin{4};
                prev_vars = varargin{5};                
            else
                error('wrong parameters')
            end
            
            if (~isscalar(val))||(~isscalar(err))
                error('only scalars are supported')
            end
            
            obj.symbolic = sym(name);
            obj.mean_val = val;
            obj.stderr = err;
            obj.units = sym(unit);
            obj.derived_from = prev_vars;
        end % c-tor
                    
        function res = char(obj)
            if sum(size(obj))>2
                res = arrayfun(@(x)(strcat([char(x) ' ; '])), obj, 'UniformOutput', false);
            else
                res = strcat([num2str(double(obj.mean_val)) '+-' num2str(double(obj.stderr)) '[' char(obj.units) ']']);
            end
        end
            
        function disp(obj)
            c = char(obj);
            if iscell(c)
                disp(['     ' c{:}])
             else
                disp(c)
            end
        end % disp
        
        function bool_res = eq(obj, val)
            if ~isa(val, 'PhysicalVal')
                bool_res = 0;
            else
                bool_res = boolean(sum(char(obj.symbolic)~=char(val.symbolic)));
            end
        end
        
%         function mean_val = get.mean_val(obj)
%             if numel(obj)~=1
%                 mean_val = arrayfun(@(x)(x.mean_val),obj);
%             else
%                 mean_val = obj.mean_val;
%             end
%         end
        
        function stderr = get.stderr(obj)
            if numel(obj)~=1
                stderr = arrayfun(@(x)(double(x.stderr)),obj);
            end
            if obj.stderr~=obj.stderr % nan
                used_vars = symvar(obj.symbolic);
                if numel(used_vars)==1
                    if isequaln(used_vars(1),obj.symbolic)
                        return
                    end
                end
                result_stderr = sym('0');
                for v=used_vars
                    partial_der = diff(obj.symbolic,v);
                    for vv=symvar(partial_der)
                        tmp = obj.derived_from(char(vv));
                        partial_der = subs(partial_der, vv, tmp.mean_val);
                    end
                    tmp = obj.derived_from(char(v));
                    %                 catch
                    %                     disp(obj)
                    %                     disp(v)
                    %                     disp(obj.symbolic)
                    %                 end
                    result_stderr = result_stderr + (partial_der*tmp.stderr)^2;
                    %                 clear tmp
                end
                obj.stderr = result_stderr^0.5;
                stderr = double(obj.stderr);
            else
                stderr = double(obj.stderr);
            end
        end
     
        function res = plus(obj1, obj2)
            if ~isa(obj2, 'PhysicalVal')
                obj2 = PhysicalVal(double(obj2), 0, obj1.units);
            end
            if ~isa(obj1, 'PhysicalVal')
                obj1 = PhysicalVal(double(obj1), 0, obj1.units);
            end
            if sum(size(obj1)~=size(obj2))
                error 'obj1 and obj2 should be the same size'
            end
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))+obj2(ind2sub(size(obj2),i));
                end
            else
                if sum(char(obj1.units)~=char(obj2.units))~=0
                    error('a and b should be the same units')
                end
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                if size(obj2.derived_from,1)==0
                    obj2_derived_from = containers.Map(char(obj2.symbolic),obj2);
                else
                    obj2_derived_from =obj2.derived_from;
                end
                res = PhysicalVal(obj1.symbolic+obj2.symbolic, obj1.mean_val+obj2.mean_val, nan, obj1.units, [obj1_derived_from; obj2_derived_from]);
            end
        end
        
        function res = minus(obj1, obj2)
            if ~isa(obj2, 'PhysicalVal')
                obj2 = PhysicalVal(double(obj2), 0, obj1.units);
            end
            if ~isa(obj1, 'PhysicalVal')
                obj1 = PhysicalVal(double(obj1), 0, obj1.units);
            end
            if sum(size(obj1)~=size(obj2))
                error 'obj1 and obj2 should be the same size'
            end
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))-obj2(ind2sub(size(obj2),i));
                end
            else
                if sum(char(obj1.units)~=char(obj2.units))~=0
                    error('a and b should be the same units')
                end
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                if size(obj2.derived_from,1)==0
                    obj2_derived_from = containers.Map(char(obj2.symbolic),obj2);
                else
                    obj2_derived_from =obj2.derived_from;
                end
                res = PhysicalVal(obj1.symbolic-obj2.symbolic, obj1.mean_val-obj2.mean_val, nan, obj1.units, [obj1_derived_from; obj2_derived_from]);
            end
        end
        
        function res = times(obj1, obj2)
            if ~isa(obj2, 'PhysicalVal')
                obj2 = PhysicalVal(double(obj2), 0, 1);
            end
            if ~isa(obj1, 'PhysicalVal')
                obj1 = PhysicalVal(double(obj1), 0, 1);
            end
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i)).*obj2(ind2sub(size(obj2),i));
                end
            else
                if sum(size(obj1)~=size(obj2))
                    error 'obj1 and obj2 should be the same size'
                end
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                if size(obj2.derived_from,1)==0
                    obj2_derived_from = containers.Map(char(obj2.symbolic),obj2);
                else
                    obj2_derived_from =obj2.derived_from;
                end
                res = PhysicalVal(obj1.symbolic*obj2.symbolic, obj1.mean_val*obj2.mean_val, nan, obj1.units*obj2.units, [obj1_derived_from; obj2_derived_from]);
            end
        end
        
        function res = mtimes(obj1, obj2)
            if isscalar(obj1) && (~isscalar(obj2))
                for i=1:numel(obj2)
                    res(ind2sub(size(obj2),i))=obj2(ind2sub(size(obj2),i))*obj1;
                end
            elseif isscalar(obj2) && (~isscalar(obj1))
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))*obj2;
                end
            elseif sum(size(obj1)~=size(obj2))
                error 'obj1 and obj2 should be the same size'
            elseif sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))*obj2(ind2sub(size(obj2),i));
                end
            else
                if ~isa(obj2, 'PhysicalVal')
                    obj2 = PhysicalVal(double(obj2), 0, 1);
                end
                if ~isa(obj1, 'PhysicalVal')
                    obj1 = PhysicalVal(double(obj1), 0, 1);
                end
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                if size(obj2.derived_from,1)==0
                    obj2_derived_from = containers.Map(char(obj2.symbolic),obj2);
                else
                    obj2_derived_from =obj2.derived_from;
                end
                res = PhysicalVal(obj1.symbolic*obj2.symbolic, obj1.mean_val*obj2.mean_val, nan, obj1.units*obj2.units, [obj1_derived_from; obj2_derived_from]);
            end
        end
        
        function res = rdivide(obj1, obj2)
            if ~isa(obj2, 'PhysicalVal')
                obj2 = PhysicalVal(double(obj2), 0, 1);
            end
            if ~isa(obj1, 'PhysicalVal')
                obj1 = PhysicalVal(double(obj1), 0, 1);
            end
            if sum(size(obj1)~=size(obj2))
                error 'obj1 and obj2 should be the same size'
            end
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i)).\obj2(ind2sub(size(obj2),i));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                if size(obj2.derived_from,1)==0
                    obj2_derived_from = containers.Map(char(obj2.symbolic),obj2);
                else
                    obj2_derived_from =obj2.derived_from;
                end
                res = PhysicalVal(obj1.symbolic/obj2.symbolic, obj1.mean_val/obj2.mean_val, nan, obj1.units/obj2.units, [obj1_derived_from; obj2_derived_from]);
            end
        end
        
        function res = mrdivide(obj1, obj2)
            if ~isa(obj2, 'PhysicalVal')
                obj2 = PhysicalVal(double(obj2), 0, 1);
            end
            if ~isa(obj1, 'PhysicalVal')
                obj1 = PhysicalVal(double(obj1), 0, 1);
            end
            if sum(size(obj1)~=size(obj2))
                error 'obj1 and obj2 should be the same size'
            end
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))\obj2(ind2sub(size(obj2),i));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                if size(obj2.derived_from,1)==0
                    obj2_derived_from = containers.Map(char(obj2.symbolic),obj2);
                else
                    obj2_derived_from =obj2.derived_from;
                end
                res = PhysicalVal(obj1.symbolic/obj2.symbolic, obj1.mean_val/obj2.mean_val, nan, obj1.units/obj2.units, [obj1_derived_from; obj2_derived_from]);
            end
        end
        
        function res = ldivide(obj1, obj2)
            res = obj2./obj1;
%             res = recalc_stderr(res);
        end
        
        function res = mldivide(obj1, obj2)
            res = obj2/obj1;
%             res = recalc_stderr(res);
        end
        
        function res = power(obj1, mantice)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))^mantice;
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(obj1.symbolic^mantice, obj1.mean_val^mantice, nan, obj1.units^mantice, obj1_derived_from);
            end
        end
        
        function res = mpower(obj1, mantice)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=obj1(ind2sub(size(obj1),i))^mantice;
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(obj1.symbolic^mantice, obj1.mean_val^mantice, nan, obj1.units^mantice, obj1_derived_from);
            end
        end
        
        function res = uminus(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=-obj1(ind2sub(size(obj1),i));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(-obj1.symbolic, -obj1.mean_val, nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res = log(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=log(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(log(obj1.symbolic), log(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res=sin(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=atan(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(sin(obj1.symbolic), sin(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res=cos(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=coc(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(cos(obj1.symbolic), cos(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res=tan(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=tan(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(tan(obj1.symbolic), tan(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res=asin(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=asin(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(asin(obj1.symbolic), asin(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res=acos(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=acos(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(acos(obj1.symbolic), acos(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
        function res=atan(obj1)
            if sum(size(obj1))~=2
                for i=1:numel(obj1)
                    res(ind2sub(size(obj1),i))=atan(obj1(ind2sub(size(obj1),i)));
                end
            else
                if size(obj1.derived_from,1)==0
                    obj1_derived_from = containers.Map(char(obj1.symbolic),obj1);
                else
                    obj1_derived_from =obj1.derived_from;
                end
                res = PhysicalVal(atan(obj1.symbolic), atan(obj1.mean_val), nan, obj1.units, obj1_derived_from);
            end
        end
        
    end % methods
end % classdef
